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Abstract 

We incorporate hydrodynamic interactions (HI) in a coarse-grained and structure-based model 
of proteins by employing the Rotne-Prager hydrodynamic tensor. We study several small proteins 
and demonstrate that HI facilitate folding. We also study HIV-1 protease and show that HI make 
the flap closing dynamics faster. The HI are found to affect time correlation functions in the 
vicinity of the native state even though they have no impact on same time characteristics of the 
structure fluctuations around the native state. 
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I. INTRODUCTION 



Globular proteins acquire distinct compact native conformations in water as a result of the 
hydrophobic effect. Another role of water is to mediate the hydrodynamic interactions (HI) 
between moving amino acids in analogy to the HI in polymers^ and between particles in 
colloidal suspensions^^ 1 ^. One expects that the HI should generally enhance cooperative 
features in the dynamics of proteins, but it is not clear in what way, exactly, will this show. 
In the case of mechanically-induced unfolding, the HI have been found to lead to a) reduction 
in peak unfolding forces when stretching at high steady velocities^, b) reduction in unfolding 
times when stretching at constant forced because of the dragging effect, and c) hindering 
of unravelling imposed through uniform^ and shear— fluid flows because of the screening 
effects. 

In this paper, we focus on the kinetics of folding and of fluctuational motions around 
the native state. We assess the relevance of HI in these phenomena within the previously 
used coarse-grained implicit-solvent model of proteins^^L^i G f N beads with the Lennard- 
Jones contact interactions. As in refs.- 1 ^, the HI are introduced through the Rotne-Prager 
tensor— 

Kikuchi et al.— have taken another approach to introduce the fluid-related effects. Instead 
of employing the tensorial field they make use of the stochastic rotation dynamics^. They 
have demonstrated that HI facilitate the collapse transition of a self-attracting homopolymer 
because of dragging effect in which a bead attracted to another bead through, say, the 
Lennard- Jones potential drags the fluid containing other beads with it. They have also, 
together with a related thesis work of Ryder—, reported a small, of order 10%, reduction 
in the folding time, tf a id, in simple models of several proteins, such as 2ci2. On the other 
hand, Baumketner and Hiwatari22 claim otherwise, pointing out that HI give rise to the 
effective repulsion between two beads which are coming towards each other, thus slowing 
down the collapse. They consider a simplified "minimal" model introduced in ref.— . For a 
short /5-hairpin, they obtain a certain increase in the folding time and lack of any effect for 
a short a-helix. 

Here, we consider four short proteins, 112y (iV=20), lbba (iV=36), lcrn (iV=46) 2ci2 
(iV=65) and one /5-hairpin (iV=14) - a fragment of the lgbl protein. In each case, our 
results are consistent with the picture of Kikuchi et al.— . The collapse phase of the folding 
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process is indeed faster. However, we find the reduction in tf ia to be more substantial than 
in the work of Ryder—, for instance, by a factor, /, of 2 for lcrn. Possible reasons for this 
quantitative discrepancy include differences in the shape of the initial conformations (Ryder 
may have more compact starting conformations compared to the fully extended situations 
considered here) and in the precise nature of the folding criterion (Ryder may have a criterion 
based on a characteristic size of a distance deviation from the native state as opposed to the 
contact-based criterion we use). Our results differ qualitatively with that of ref.—. It should 
noted that the model used there comes with attractive contacts only between hydrophobic 
amino acids and requires this attraction to be the strongest at the distance of 4.2 A between 
the corresponding C a 's. In proteins, contact distances between the C a atoms extend in fact 
from 4.4 to 12.8 A— which is likely to change the nature of the effects associated with the 
HI. In fact, when we a) constrain the potentials for the freto-hairpin case to have a minimum 
at 4.2 A the reduction factor goes down to 1.2 and then, in addition, b) increase the number 
of allowed contacts to imitate the situation envisioned within the minimal model, / becomes 
equal to 0.87, indicating that the HI make the folding longer. Thus the result of ref.— is 
related to the minimal nature of the model considered there. 

We then consider a larger protein - the HIV-1 protease (la30). This protein is a ho- 
modimer and each of the two chains comprises 98 amino acids. Its active site is covered 
by two flexible /3-hairpins, called flaps, that controll the entry of a polypeptide substrate. 
The flap opening dynamics in this protein has been studied before within all atorn^i and 
coarse-grained molecular dynamics schemes. Here, we pull the flaps apart and then monitor 
their return to the native form as a function of time. We observe that the HI make the flap 
closing faster. 

Finally, we return to the small protein 112y and investigate fluctuations around the native 
state by considering single- and double-residue characteristics of these fluctuations. We find 
that the HI have no impact on their same-time averages. This result is not surprising because, 
in the overdamped limit considered here, the equilibrium distributions of conformations are 
the same. However, it sets the stage for the observation that time correlation functions are 
sensitive to the HI despite the protein being essentialy folded. 
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II. METHODS 



The coarse-grained, Go-type model 2 ^ of a protein we use is constructed based on the 
knowledge of the native state. There are many ways to implement it. Examples are given in 
refs.— i22j22j2Ij22j22. The details of the particular implementation we use are described ini 4 -^. 
Each residue is represented by a single bead centered on the position of the C a atom. The 
successive beads along the backbone are tethered by harmonic potentials with a minimum at 
3.8 A and they are also endowed with the chirality based local backbone stiffness 14 . The other 
interactions between the residues i and j are split into two classes: native and non-native 
as described in ref.— . The native contacts are endowed with the effective Lennard- Jones 



potential Vu = 4e 
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whereas the non-native contacts are purely repulsive. 
The length parameters are chosen so that the potential minima correspond, pair-by-pair, 
to the experimentally established native distances between the respective aminoacids. The 
energy parameter, e, is taken to be uniform and 0.3e/fce is usually playing the role of the 
room temperature^. 

The time evolution is described by the Brownian dynamics (BD) 34 . The displacement of 
particle i in time step At obeys 

- r? = £(V, • D° ) At + ^ E D £ ' F ? A * + B - W 

where the index denotes the values of respective quantities at the beginning of the time 
step, Fj is the force exerted on particle j by other particles, and T is temperature. D is a 
diffusion tensor and B - a random displacement given by a Gaussian distribution with an 
average value of zero and covariance obeying 



< BiB,- >= 2D°-At. 



(2) 



If the diffusion tensor is nondiagonal, there exists a hydrodynamic coupling between particles 
% and j (cf. EqJT]). We take the diffusion tensor in the form 1 ^ 1 ^ 
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where = rj — r^, a is the hydrodynamic radius, and 7 is the damping constant. By 
Stokes's law, 7 = 6irarj, where rj is the viscosity of the fluid. The characteristic time scale in 
the problem, r, is of order 1 ns, which reflects duration of a diffusional passage of a typical 
contact distance (~ 5 A). 

The situation corresponding to the lack of HI will be denoted by BD (for Brownian 
dynamics) and will be implemented by setting (i 7^ j) to zero. The selection of the 
value of the hydrodynamic radius is not obvious and more studies are needed to settle this 
issue. We just want to determine qualitatively what would happen if it was non-zero. When 
one thinks of amino acids as represented by spheres then one would expect that a should not 
exceed a half of the distance between consecutive C a 's. Thus we consider the hydrodynamic 
radius of 1.5 A to have a substantial enough and yet satisfying this criterion. On the other 
hand, the side groups extend laterally and may produce a correspondingly larger drag force. 
It has been argued 3 ^ 1 ^ that a characteristic a of an amino acid could be even as big as 4.2 A 
whereas a characteristic van der Waals radius is about 3.0 A—. The van der Waals volume 
does not include hydratation layers. Such big values would mean existence of an overlap 
between spheres in a chain, but its usage takes into account the non-spherical properties of 
the amino acids in an approximate way. Considering a of 1.5 A has a numerical advantage 
because the time unit is governed by the damping constant and is thus proportional to 
a. Thus the folding times are expected to scale linearly with a for systems both with and 
without HI, making it relevant to assess the role of HI regardless the particular choice of 
a. It cannot be ruled out, however, that HI may introduce some corrections to scaling. 
Assessment of such corrections would need a separate study. 

It should be noted that there are some amino acid to amino acid variations in the values 
of the van der Waals and the de la Torre - Bloomfield 3 ^ hydrodynamic radii. For glycine, 
alanine, and arginine, the former are 2.4 A, 2.6 A, and 3.3 A respectively, and the latter 
are 3.6 A, 3.7 A, and 4.5 A respectively 3 ^. These variations are not expected to affect the 
findings in analogy to the lack of sensitivity to the variations in the values of the masses^ 3 - 
or in the values of the damping constant 39 . A scaling by the mean value of the parameters 
has been demonstrated in these other cases. 

Throughout this paper we employ the Brownian dynamics evolution algorithm since it 
allows for a straightforward incorporation of the HI. However, all of the BD results obtained 
here are reproduced by the Langevin dynamics algorithm as used in refs.— The Langevin 
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dynamics involves inertia but the system is overdamped. 

III. RESULTS 

A. Kinetics of folding 

Figure 1 shows the T-dependence of the median value of tf Q u determined as the first time 
to establish all native contacts (r^ < 1.5ery). For each of the proteins studied, there is a 
broad plateau of optimal folding. In the optimal regime, tf id with HI is shorter by a factor 
/ than in the BD case. The value of the factor depends on the protein and is indicated on 
the right hand side of the figure. For the /3-hairpin, / is 1.3. The largest protein studied, 
2ci2, has been simulated only at one temperature (0.3 e/k B ). The simulations yielded tf u of 
~ 840r and ~ 470r for the BD and HI models respectively. The corresponding / factor of 1.8 
is clearly distinct from ~1.1 found in ref.—. The discrepancy may probably be attributed 
to two circumstances: a) the statistics involved in ref.— have been restricted to several 
trajectories, b) the values of / depend on the nature of the starting conformations. The 
bigger the number of the native contacts that are present in the initial stages of folding, the 
smaller the value of /. In our simulations, the starting conformations are fully extended. 

The Hi-induced acceleration of folding is governed mostly by the initial collapse as il- 
lustrated for two typical trajectories for lcrn in the bottom panel of Figure 2. This figure 
shows the time evolution of the radius of gyration, R g . Qualitatively, the collapsed phase 
is said to be reached when R g does not exceed 15% of its native value. In agreement with 
ref.—, the HI are seen to significantly accelerate arrival at the collapsed phase (by a factor 
larger, ~ 2.6, than the acceleration factor found for t/ ^). The subsequent establishment of 
contacts involves a stochastic search in the conformational space and is not affected by any 
hydrodynamic effects. 

The top panel of Figure 2 shows the distributions of the values of tf o[d across 100 trajec- 
tories for lcrn. The distributions are non-Gaussian and have pronounced tails. The peak 
of the distribution corresponding to HI is at shorter times than BD, but the long-time ends 
are comparable. It should be noted that both at high and low temperature ends, when the 
folding time becomes large, the distinction between the HI and BD cases evaporates because 
when the protein kinetics are slow the related induced fluid flows are sluggish. 
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Figure 3 shows the scenario diagram 14 which represents the first times needed to establish 
specific native contacts as averaged over the trajectories. The native contacts are labelled by 
their sequential separation \j — i\. It is seen that the presence of HI accelerates establishment 
of all contacts without changing the order in which they are set. The same statement applies 
to the other proteins we have studied and we expect it to hold in general. 

B. Closing of flaps in HIV-1 protease 

The native structure of HIV-1 protease is schematically represented by the lower right 
conformation shown in Figure 4. The hairpin fragments that form the flaps are highlighted 
by showing the corresponding C a atoms. The locations of the flaps can be conveniently 
described by providing coordinates of the four-atom centers of mass. The two centers of 
mass are _R/ jnat =4.09 A away from each other. We pull the model flaps apart by exerting a 
stretching force applied along the direction that links the initial centers of mass of the two 
flaps until they are separated by 32.75 A. The resulting conformation is shown on the upper 
left in Figure 4. The force is then removed and we monitor the distance, Rf, between the 
centers of mass as a function of time. 

Figure 4 demonstrates that the HI make the closing of the flaps faster. The value of Rf = 
2Rf nat is reached about twice as fast with HI compared to the BD case. The phenomenon 
is analogous to refolding discussed previously. However, it affects only a part of the full 
protein. 

C. Structure fluctuations around the native state 

We now consider 112y which is the smallest among the set of proteins studied here so that 
an appropriate conformational statistics can be generated. Its native structure is shown in 
Figure 5. It is akin to a /3-hairpin in wich one of the "arms", however, is shaped into an 
a-helix. The 20 amino acids are enumerated from the helical end. We set the protein in its 
native conformation, then evolve it at T=0.3e/fce in four different trajectories which last for 
200 000 r. 
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1. Equal time characteristics 



One way to quantify the dynamics of the motion around the native state is to study the 
fluctuations of amino acid positions < AR 2 (t) > around the average value. We follow the 
procedure of Kabsch and Sanders^ and project any evolved conformation onto a reference 
conformation so that the results are not affected by translation of the center of mass and 
by rotation of the whole molecule. < AR 2 > is defined then as < ff > — < r*j > 2 . The 
top panel of Figure 6 shows the modulations of < AR 2 > along the sequence. The largest 
fluctuations take place at the termini and at i=12 where the helix joins the other arm of the 
hairpin. Smaller, but still substantial, fluctuations occur at i= 8 and 15. The data points 
corresponding to the HI and BD cases nearly overlap indicating the lack of relevance of the 
HI in the equilibrium fluctuations. 

A similar statement applies to the two-point (equal time) characteristics of the motion 
as illustrated in Figure 7 for pairs of residues i and j that make a native contact. There are 
27 such contacts in 112y (counted by the index I) and 13 of them are indicated by the lines 
in Figure 5. The allocation of the index I to its corresponding pair of i and j is also written 
underneath the data point in the top panel of Figure 7. 

The simplest characteristic is fi = (< r 2 j > — < > 2 ) 1 ' 2 . This quantity does not 
depend on any rigid body motion of the protein and it rapidly acquires lack of dependence 
on the duration of the measurement. The top panel of Figure 7 shows that it is insensitive 
to the inclusion of the HI. The modulations in the strenghts of fi span a factor of two. Large 
values indicate high amplitude oscillations (like in contact 22 between amino acids 7 and 
11) and small values point to a relative rigidness (like in contact 27 between 11 and 14). A 
removal of the "keystone" contact 18 (between residues 6 and 17) is found to reduce ^22 and 
/15 but to unhance /s and f$. For an isolated a-helix the f]s are usually small and uniform 
except for enhancements at the termini. The /3-hairpin has largest fluctuations at contacts 
that link the termini (fi ~ 1.5 A). 

Another popular characteristic is the so called dynamical cross-correlation map (see e.g. 
refs.— ^i^ 2 -) defined as 

Q = C tj = <Ar t - Ar, > /[< Ar 2 > < Ay) >] 1/2 . (5) 

Unlike fi, C\ involves features which are primarily orientational and does not depend much 
on the distance between i and j. C\ can be either positive or negative, depending on 
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the nature of the relative displacement. Its determination requires making the rigid-body 
transformation that brings the conformation closest to the reference structure. We find that 
C\ for U2y settles in their stationary values within 2000 r. These stationary values are 
shown in the bottom panel of Figure 7. The values of Ci are modulated but the sensitivity 
to the presence of HI is again seen to be minor. 

All of the above single- and double-residue equal time characteristics correlation functions 
may be also calculated by using the Gaussian network model^ and by performing the normal 
mode analysis. We find that the nature of the sequence- and pair-dependent modulations 
coming from this approach is as in Figures 6 and 7 respectively. 

2. Time- dependent characteristics 

We now take various conformations along each trajectory and investigate what happens 
to them a time t later. The motion of the center of mass of the protein is subtracted and 
we could get reasonable statistics for time intervals not exceeding 5000 r. We define the 
correlation function 

Si(t) =< fi(0) ■ n(t) > I < r-(0) ■ f,(0) > (6) 

and S(t) =< Si > av , where < ... > denotes an average over the trajectories and < ... > av 
an average over the residues. The definition implies that for t=0, Si—1. The initial decay 
of S(t) with t is shown in Figure 8. It is seen that the decay, or a 'decoherence', from the 
initial state is faster, when the HI are included. Due to limitations in statistics at longer 
delay times, we cannot pinpoint the precise functional law for the time dependence of S(t). 
However, the important observations is that the time dependence of time correlations is 
affected by the presence of the HI. 

It should be noted here that there exists a puzzle regarding the origin of long-time tails 
in the correlation functions in proteins. Namely, the experiments by Xie and coworkers^ 
seem to show that fluctuations in proteins decay very slowly in time, following a power law. 
A number of theoretical attempts to explain this time-dependence have been made^^^ 1 ^ 
but most of them predict the decay on a much shorter time-scales than those measured in 
experiments, or perhaps they do not take into account the presence of cross links in the 
protein. 
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A related time-dependent characteristic is provided by 

Ai(t) =< \m - rl(0)] 2 > / « \n{t) - rl(0)] 2 » av (7) 

in which the numerator would be related to the single-residue diffusion coefficient, Di(t), if 
divided by 6t. However, this D^t) decays to zero, instead of homing in on a constant, since 
the motion of the center of mass is excluded in the trajectory. By dividing the numerator 
by the denominator in the expression for A{(t), one effectively removes most of the time 
dependence, whether it is provided by BD or by HI. The interesting observation is, as 
demonstrated in the bottom panel of Figure 6, that Aj(£) shows no difference between BD 
and HI and, practically, has no time dependence. Its sequential behavior is qualitatively 
similar to that of < ARf > (the top panel of Figure 6). 

In summary, the HI affect the time scale of folding significantly by making the collapse 
faster through the cooperative action of the drag forces. This phenomenon can be equiv- 
alently described by invoking a system with no HI but with a reduced effective viscosity 
by about a factor of two. Thus in the context of protein folding, or phenomena similar to 
the flap closing in HIV-1 protease, the presence of HI just shifts the effective value of the 
damping constant 7. Observing phenomena that can be clearly attributed to the HI may 
then be difficult experimentally unless one investigates a broad range of temperatures. In 
an equilibrium evolution, the HI affect time correlation functions without influencing equal 
time correlations such as the rms fluctuations in the contact length. 

This paper has benefited from many fruitful discussions with P. Szymczak who also 
provided us with his numerical code. We also appreciate discussions with A. Giansanti and 
M. Wojciechowski. This work has been supported by the grant N N202 0852 33 from the 
Ministry of Science and Higher Education in Poland. 
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FIGURE CAPTIONS 



Fig. 1. The median folding times as a function of temperature for the three proteins. The 
dotted lines and open data points correspond to the BD model whereas the solid lines 
and data ponts to the model with the HI. The hydrodynamic radius is 1.5 A. The 
/-factors on the right give the ratios of the optimal folding times without the HI to 
those with the HI. 

Fig. 2. Top panel: The distribution of single trajectory folding times at k B T /e=0.3 for 
the model crambin. Bottom panel: Examples of time evolution of the radius of gyra- 
tion. The down-pointing arrows indicate entries into the phase of compact collapsed 
conformations. The up-pointing arrows show arrivals at the native conformation. 

Fig. 3. The scenario diagram for folding of the model lcrn. The hexagons correspond to 
the time to establish the constacts between C and I for the first time. The circles 
correspond to the contacts between two /3-strands. The full circles - to the contacts 
within helices and between helices. The asterisks correspond to all other contacts. 
The upper symbols involve the HI and the lower do not. 

Fig. 4 The flap dynamics in la30 as discussed in the main text. The conformation on the 
right is native and on the left is stretched. The dotted line corresponds to the BD 
case and is based on 50 trajectories. The solid line corresponds to the HI case and is 
based on 20 trajectories, all starting from the same stretched state. 

Fig. 5. The native state of 112y. The larger numerical symbols enumerate the amino acids 
from the N terminal. The smaller symbols enumerate native contacts (the C-terminal 
generates no native contacts). The contacts are also indicated by lines. The solid 
lines correspond to contacts with strong distance fluctuations, as shown in Figure 6, 
whereas the dashed lines to those with weak fluctuations. 

Fig. 6. Top panel: Variance in single-residue position fluctuations as enumerated se- 
quentially. The solid (open) symbols are for the model with (without) the HI and 
/c#T/e=0.3. Bottom panel: normalized time-dependent variance Aj(t) for t = lOOOr. 
The symbols corresponding to HI and BD in this and following two figures are dis- 
placed laterally around the proper position to enhance their individual visibility. 
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Fig. 7. Fluctuations in the native contacts around the native state at ksT /e=0.3. Top 
panel: the distance rms fluctuations. The numbers indicate pairs of amino acids that 
are are connected by the native contact labeled by /. Bottom panel: the orientational 
rms fluctuations (known also as the dynamical cross-correlations). 

Fig. 8. The normalized time correlation S(t) as a function of time. 
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